﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "ingamma")]
    public static unsafe double ingamma_dll(double a, double x)
    {
        return ingamma(a, x);
    }

    /// <summary>
    /// 不完全Gamma函数
    /// </summary>
    /// <param name="a">a自变量值。要求a>0。</param>
    /// <param name="x">x自变量值。要求x>=0。</param>
    /// <returns>函数返回不完全Gamma函数值。</returns>
    public static unsafe double ingamma(double a, double x)
    {
        int n, flag;
        double p, q, d, s, s1, p0, q0, p1, q1, qq;

        if ((a <= 0.0) || (x < 0.0))
        {
            return (-1.0);
        }
        if (x + 1.0 == 1.0) return (0.0);
        if (x > 1.0e+35) return (1.0);
        q = Math.Log(x);
        q = a * q;
        qq = Math.Exp(q);
        if (x < 1.0 + a)
        {
            p = a;
            d = 1.0 / a;
            s = d; n = 0;
            do
            {
                n = n + 1;
                p = p + 1.0;
                d = d * x / p;
                s = s + d;
                flag = (Math.Abs(d) >= Math.Abs(s) * 1.0e-07) ? 1 : 0;
            } while ((n <= 100) && (flag > 0));

            //if(!flag)
            if (flag == 0)
            {
                s = s * Math.Exp(-x) * qq / gamma(a);
                return (s);
            }
        }
        else
        {
            s = 1.0 / x;
            p0 = 0.0;
            p1 = 1.0;
            q0 = 1.0;
            q1 = x;
            for (n = 1; n <= 100; n++)
            {
                p0 = p1 + (n - a) * p0;
                q0 = q1 + (n - a) * q0;
                p = x * p0 + n * p1;
                q = x * q0 + n * q1;
                if (Math.Abs(q) + 1.0 != 1.0)
                {
                    s1 = p / q; p1 = p; q1 = q;
                    if (Math.Abs((s1 - s) / s1) < 1.0e-07)
                    {
                        s = s1 * Math.Exp(-x) * qq / gamma(a);
                        return (1.0 - s);
                    }
                    s = s1;
                }
                p1 = p; q1 = q;
            }
        }
        s = 1.0 - s * Math.Exp(-x) * qq / gamma(a);
        return (s);
    }

    /*
    // 不完全Gamma函数例
      int main()
      { 
          int i,j;
          double y,s,t;
          double a[3]={0.5,5.0,50.0};
          double x[3]={0.1,1.0,10.0};
          for (i=0; i<=2; i++)
          for (j=0; j<=2; j++)
          { 
              s=a[i]; t=x[j];
              y = ingamma(s,t);
              cout <<"ingamma(" <<a[i] <<", " <<x[j] <<") = " <<y <<endl;
          }
          return 0;
      }
    */
}

